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FOREWORD 


This is a progress report on the research project "Numerical Solutions 
of Three-Dimensional Navier-Stokes Equations for Closed-Bluff Bodies," for 
the period ending June 30, 1987. Specific attention has been directed to 
investigate the "Conservative Finite Volume Solutions of a Linear Hyperbolic 
Transport Equation in Two and Three Dimensions Using Multiple Grids." This 
work was supported by the NASA/Langley Research Center (Computer Appli- 
cations Branch, Analysis and Computation Division) through the cooperative 
agreement NCC1-68. The cooperative agreement was monitored by Dr. Robert E. 
Smith of the Computer Applications Branch, ACD. 
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CONSERVATIVE FINITE VOLUME SOLUTIONS OF A LINEAR HYPERBOLIC TRANSPORT 
EQUATION IN TWO AND THREE DIMENSIONS USING MULTIPLE GRIDS 


By 

M. Kathong 1 and S.N. Tiwari 2 

ABSTRACT 

The feasibility of the multiple grid technique is investigated by 
solving linear hyperbolic equations for simple two- and three-dimensional 
cases. The results are compared with exact solutions and solutions obtained 
from the single grid calculations. It is demonstrated that the technique 
works reasonably well when two grid systems contain grid cells of compar- 
ative sizes. The study indicates that use of the multiple grid does not 
introduce any significant error and that it can be used to attack more com- 
plex problems. 


Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23529. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23529. 
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1. INTRODUCTION 


In recent years much progress has been made in the solutions of steady 
state equations of motion in both two and three dimensions. For complex 
shapes these calculations are usually based on a body-fitted curvilinear 
grid. For general three-dimensional bodies (for example, an aircraft 
configuration), it is very difficult to construct a body-fitted coordinate 
system [1]*. To simplify this problem, it is becoming more common to use 
several grids (multiple grids) at once, each in a different coordinate 
system [2-4]. This approach results in new boundaries within the given 
region at the interfaces of the various grids (Fig. 1.1). In order that 
information be transferred from one grid to another accurately, it is impor- 
tant to treat grid points on the interfaces with care. The non-linear 
nature of the equations of motion permits solutions with discontinuities 
such as shocks and slip surfaces. In order that such discontinuities assume 
the right strength and physical location, it is imperative that the scheme 
used for the calculations be conservative [5]. In a multiple grid calcu- 
lation, it is important that the interfaces are also treated in a conser- 
vative manner so that the discontinuities can move freely across these 
i nterf aces. 

The question of conservation when switching between two different grids 
or numerical schemes has been considered by several authors. Warming and 
Beam [6] derived transition operators for switching conservatively between 
MacCormack's method and a second order upwind scheme. Hessenius and Pulliam 
[7] applied this transition operator approach to derive so-called zonal 


*The nimbers in brackets indicate references. 
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interface conditions, but with a significant loss of accuracy at the zonal 
interfaces. Rai [8] has developed conservative zonal interface conditions 
for zonal grids which share a common grid line, and has nice calculations 
demonstrating the shock capturing ability of zonal grids with a discontinu- 
ity crossing zones. Earlier work in the area of multiple grids includes 
that of Cambier et _al_. [9] who analyzed the zonal -boundary problem for a 
system of hyperbolic equations and used the compatibility equations to de- 
velop a zonal-boundary scheme. Good results are presented for transonic 
channel flow. However, the use of the compatibility equations results in a 
zonal-boundary scheme that is not conservative and, hence, unsuitable for 
problem in which flow discontinuities move from one grid to another. Rai 
et _al_. [10] present results obtained on metric discontinuous grids; the 
integration scheme used is the Osher upwind scheme. Reference 11 gives 
results obtained on overlaid grids in conjunction with the stream function 
approach. For more information on multiple grid one should refer to [12- 
20]. the need for conservative grid interfaces is illustrated by Benek, 
Steger, and Doughtery [3]. 

Dukowicz [21] described a rezoning method for arbitrary quadrilateral 
meshes in two dimensions. Ramshaw [22] suggested a procedure which is sim- 
ilar to the method of Dukowicz, but is simpler and more direct. Note that 
they defined rezoning as a method for transferring a conserved quantity, Q, 
from one generalized mesh to another when the volumetric density of Q is 
uniform within each cell of the original mesh. A computer program which 
follows Ramshaw’ s procedure has been written and tested with various types 
of grids and variables. The program seems to be working well. The objec- 
tive of this study is to establish whether or not this technique is feasible 
for various grid systems while keeping the governing equation simple. 
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For simplicity, the scalar 2D and 3D hyperbolic equations are chosen 
for this study. Hyperbolic equations have had a history of being used as 
model equations for testing newly developed schemes, for example, see [23- 
30], Finite volume approach is chosen along with three-stage Runge-Kutta 
time integration. The three-stage Runge-Kutta integration is of 2nd order 
accuracy. Many authors have applied the finite volume approach and Runge- 
Kutta to their numerical calculations [38-42]. 

This study can be divided into two portions. First, the scalar 2D 

hyperbolic equations, q t + aq + bq = 0 is considered. Here, the equation 

t x y 

is solved on a two dimensional grid system which is changed into another 
grid system at some time, say t = t^. The information obtained from the 
first grid calculation is transferred to the second grid by Ramshaw's tech- 
nique. Many authors, for example Berger [43], suggested that in order to 
retain the conservative property of the numerical scheme, the interpolation 
of flux across the interface is needed. For this problem, q itself is the 
flux. Ramshaw's technique is used to transfer q across the interface. 

The second portion of this study is to examine the technique with the scalar 
3D hyperbolic equation, q f + aq + bq + cq_ = 0. This equation is more 
suitable as the model equation of the equations of motion. Here, the grid 
system is changed from one to another at some x = y plane, say plane z = z^. 
Plane z^ is the interface mentioned previously. For this problem, the flux 
is going through this plane, i.e., flux in z direction = h = cq, needs to be 
transferred across the plane. Again, Ramshaw's technique is implemented in 
doing so. Details and procedures of solving both 2D and 3D hyperbolic 
equation is given in later sections. The results are compared with exact 
solutions and solutions obtained from the single grid calculations, i.e., 
without applying the Ramshaw's technique. It is important to note that 


4 



the objective of this study is to determine whether the technique is appli- 
cable for these simple governing equations. If it is applicable, one may 
have some confidence that it will be applicable for more complicated equa- 
tions. It is expected that the results from this study will yeild a signi- 
ficant contribution in the area of Computational Fluid Dynamics. 

In Sec. 2, a brief discussion on Ramshaw's technique is given. Sec- 
tions 3 and 4 are concerned about the two and three dimensional hyperbolic 
equation, respectively. Finally, the conclusion is given in Sec. 5. 

2. DISCUSSION ON CONSERVATIVE REZONING ALGORITHM 

Rezoning is a method for transferring a conserved quantity Q from one 
generalized mesh to another when the volumetric density of Q is uniform 
within each cell of the orginial mesh. The need for rezoning was discussed 
in Sec. 1. A rezoning method for arbitrary quadrilateral meshes in two 
dimensions has recently been described by Dukowicz [21]. Ramshaw [22] 
suggested the procedure in doing so which is similar in spirit to the method 
of Dukowicz, but is simpler and more direct. A computer program following 
this procedure has been written and is working well for example grids and a 
wide variety of choice of variables. 

By far the most common type of generalized mesh is the arbitrary quad- 
rilateral mesh, which is convenient to work with because it has the same 
simple topological and logical structure as a square or rectangular mesh, 
the basic idea behind the rezoning is simple. Here, two grid systems are 
overlapped each other in some fashion (Fig. 2.1). The conserved quantity, 
denoted by Qo^, is to be transferred from the old grid system (Ao^ is the 
area of each mesh in the old grid system) to the new grid system in which 
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An,, is the area of each. The quantity Qn. , is denoted as the transferred 

J 1 J 

quantity in each mesh of the new grid system. Thus, Qn,, can be computed 

1 J 

by: 


Q n i j = E(Qo kl ) (Ano k1 ) = (qo k] ) (Ano k] ) 
Ao ki 


where An„ is the portion of the area An., which is contained in the area 
o r ’ 1 


kl ^ 

Ao^ , and the summation is up to the number of the old meshes contained in 

An.j. The quantity q^ = Qo kl represents the volumetric density of Qo kl 

and is assumed to be uniform within each cell. The task, now, is to find 

Ano, , and the number of the old meshes contained in each An. .. The area of 
kl ij 

the polygon P is given by [44] 


Aq = 1 | E (^i ~ *2 ^1 ^ * 

H J s s 

where the summation is over all the sides of P, and is either +1 or -1 

depending on whether P lies to the left or right of side s. Note that the 

s s s s 

endpoint coordinates (x^, y^) and (x^, y 2 ) are considered to be associated 
with the side s and not with the particular polygon P. 

Figure 2.2 shows arbitrary overlap grids. The overlap areas are poly- 
gons whose sides are segments of the old-mesh lines. The number of sides of 
each type, and total number of sides, will be different for different over- 
lap areas. Each side is common to two overlap areas, the one on the left 
(L) and the one on the right (R), and these overlap areas may be considered 
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to be associated with the side. 

The objective is to apportion a conserved quantity Q, whose volumetric 

density q is considered uniform within each cell of the old mesh, into the 

cells of the new mesh. It would be inefficient and difficult to automate in 

a computer, to naively sweep over the overlap areas directly. Instead, 

Ramshaw suggests to evaluate the same contributions by sweeping over the 

sides or segments s. The side or segment s is any side or segment of the 

polygon (overlap area). The coordinate of the two endpoints of side s 

s s s s 

will be denoted by (Xp y^) and (x 2 , y 2 ). 

If the side s is a segment of the old mesh, then the quantity Q in 
the new-mesh cell containing side s is to be incremented by an amount: 

A° = i (q L - q R ) (xp y 2 ) and (x|, yj) 

If the side s is a segment of the new mesh, then the contribution to cell 
on its left is A^| = — q Q (q L - q R ) (x*, y*) and (x*, y*) while the 

N 

contribution to the cell on its right is just -A s , where q Q is the 

volumetric density of the quantity Q of the old mesh cell in which side s 
o N 

lies. Adding a $ and a $ for each of the new mesh cell yields the quantity Q 
contained in each of the new mesh cells. 

3. TWO-DIMENSIONAL HYPERBOLIC EQUATION 

Consider the two-dimensional hyperbolic equation, 

q + aq + bq = 0, 0 <x<l, 0<y<l (3.1a) 

t x y 
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with the initial condition. 


q(x,y,0) = f(x,y) (3.1b) 

The exact solution of Eq. 3.1 can be found as 

q(x,y,t) = f(x - at, y - bt) (3.2) 

This problem is well posed [45] if boundary conditions are specified at 
the boundaries where the dot product, (a,b)» n > 0; here (a,b) = a"e + 

A 

b£ and 'n' is the unit normal at the boundaries. This problem is solved on a 

y 

multiple grid system using the finite-volume approach along with the three- 
stage Runge-Kutta integration in time. Section 3.1 is concerned with the 
formulation of the 2D hyperbolic equation on the curvilinear grid system. 

The implementation on multiple grid system is briefly discussed in Sec. 3.2. 
Finally, the results and discussion on various grid configuration and ini- 
tial condition are given in Sec. 3.3. 


3.1 Formulation of 2D Hyperbolic Equation on Curvilinear Grid 
Recall Eq. 3.1, 

q t + aq x ♦ bq y =■ 0 


Since a and b are constants, one can write Eq. 3.1 in the conservation 
for as [46]. 

q t + f x + g y = 0 ( 3 - 2 ) 

where 


f = aq, g = bq 
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Equation 3.2 could be solved on any (x,y) coordinates, called physical 
coordinates. Instead, one normally solve Eq. 3.2 on the orthogonal 
coordinates system, called computational coordinates [47], Thus, equation 
3.2 needs to be transformed to the computational coordinates. Even though 
the transformation into the computational coordinates produces some 
additional terms in the transformed partial differential equation (thus 
increases the amount of computation required), it has some advantages. 
Boundary conditions can be treated in the straight forward way. Also, 
transformed equations are solved in the simple orthogonal and possibly equal 
spaced region. Severe departure from orthogonality will introduce 
truncation error in difference expressions [47]. Figure 3.1 illustrates the 
physical versus computational domain in a two-dimensional space. By using 
the chain rule, Eq. (3.2) is written as: 

q t + * <V* + vV *° <3 - 3) 

Also using the chain rule, one can write: 

dx = x„d? + x dn 
£ n 

dy = y r d£ + y dn 

£ n 

which can be written in the matrix form as: 


— - 





dx 

= 

x £ x n 


d? 

dy 


h y n 


d" 

_ „ 


— — 


— — 
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Physical Domain 

Fig. 3.1. Physical vs. Computational Domain 



Computational Domain 
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The inversion of the matrix yields: 


— — 1 


r— __ 

— — 


= 

x„ X 
l n 


dx 

dn 


y c y 


dy 



S n 



i— — 


M - -- 



Simi 1 arly, 


d£ = £ dx + l dy 
x y 

dn = n dx + n dy 
x y 




From Eqs . 3.4 and 3.5, it can be seen that: 



where 


Vn 


Vn 


is the Jacobian of transformation. Obviously, Eq. 3.6 exists when 



x„ x 
y £ y n 


* 0 . 


From Eq. 3.6 it is found that C x = Jy^ 5 = 


n 


x 




Jx 


K 


(3.4) 

(3.5) 

(3.6) 
(3.6a) 

J * » or 
(3.6b) 
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Equation 3.3 can be written as: 


I q t + 

1 <ft„ - 95 y > 

+ 

1 (fn + gn ) 

— x y 

J 

J 

5 

J 


q t + f 5 + 9n = 0 


(3.7) 


q = q/J 

f = I (f? x + gc y ) = fy n - gx^ (3.7a) 

J 

g = I ( fri x + g^y) = - f y c + g* c 

j 

Equation 3.7 is the transformed equation which is solved on the 
computational coordinates. 

3.1.1 Difference Equation 

In this section, the finite volume approach is applied to Eq. 3.7 along 
with three-stage Runge-Kutta integration in time. The formulation of the 
difference equations on the equally spaced in each direction (c - n), 
computational domain is given below. 

Consider a computational cell illustrated in Fig. 3.2. Recall Eq. 3.7 

q t + f 5 + q „ = 0 
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Integrating over the cell yields 


f E q t dE • - ,f E (f 5 + 9„)dE. 


Using the divergence theorem, one obtains 


/ E q t dE » -(/ c 'f { n*- dt * ■ ds) 


Applying this to Fig. 3.2, one can write 


a? An^ t 


= -An(f - f , J 

i + 1 , j + 1 i + 1, j + 1 ij + 1 


~ A ^ 9 i + i ( j +1 " g i + ^ 

2 2 


I«t 

J 


✓N /\ 

i + I , J + I = (f i + 1* J + 1 ‘ V J + I 
2 2 AC 2 2 


( " 9 i + I , j + 1 " 9 i + I , j) 

An ~ o 


<*t 


i + i » J + I 
2 2 


i + I, j + I 
2 2 


i + 1, j + i " f 1, j + 

A? 0 0 


^ + 1, 3 


+ 1 


/\ 

g 


An 


1 + I, J 
2 


t) 


f(q n ) 


(3.8) 
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Each term on the right hand side of Eq. 3.8 can be written as 


Vi*i“ (fy n - g \> i ♦ i, j ♦ i 


a ^i’ j + 1 ~ y 1,3^"^ x i* j + 1 x i’ ^ 9 1 * 3 + 1. 


'l + 1., j ^' fy 5 + 9 C* 1 ♦ J., 3 


‘ " a * y i + 1, 3 ’ y >. J* +b(x i + 1, 3 " X ’. J i *iJ 

AC 2 

q i , j + JL — ^i + 1, j + 1. + q i - l, j + 1^ q i + i» 0 

2 2 2 2 2 2 2 


" I (q i + 1, j + 1 + q i + 1, j - 1 } 
2 2 2 2 2 


X C i + 1. .i + 1 — ^ X i + 1, 3 + 1 *i» j + 1 + 1 UA 

2 2 2 AC AC 


i + 1, j + 1 = - (y i + 1, j + 1 ~ y i, j + 1 + y i + 1, j ~ y i, J ) 
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i + 1, j + 1 — ^ x i + 1, j + 1 ~ x i + 1, j + x i, j + 1 " x i, p 


2 2 2 


An 


An 


i + I* J + I — + 1, j + 1 ~ y i + 1, j + y i , j + 1 ~ y i , j ^ 

2 2 2 An An 


The exact boundary conditions, q(Xp yp t) = f(x^ - at, y^ - bt), are 
given where the product (a, b) • n > 0, and the extrapolated values of q 
are given at the boundaries where (a, b) • n < 0. This is due to the nature 
of the hyperbolic equations. 


3.1.2 Three-Stage Runge-Kutta Integration 
Recall that Eq. 3.8 is written as 


where 


q t I i + 1 , J + I 
2 2 


- W) 


q = q(t n> 1 + 1, j + 1 
2 2 


By applying three-stage Runge-Kutta, one may express: 


n + 1* n 

q 3 q 


+ Atf(q n ) 


n + 1** n 

q = q 


+ _1 Atf(q n ) + _1 Atf(q n + * ) 
2 2 


q" + 1 



+ ^Atf(q n ) +_1 Atf(q n + * ) 

2 2 
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3.2 Implementation on Multiple Grid System 


As mentioned previously, for multiple grid system, the interpolation of 
flux will yield a conservative scheme. For the scalar 2D hyperbolic 
equation, Eq. 3.1, q itself is the flux transporting in time. For this 
problem, the difference equations as described in Sec. 3.1.1 are solved on 
one grid system from initial time t = 0 up to some time, say t = t^. Then, 
at t = t^ the grid system is changed and the difference equations are solved 
on this new grid system. The initial conditions for the new grid system are 
obtained by the interpolation of the final values in the old grid system, 
i.e., q(x,y,t^). Ramshaw's technique, as described in Sec. 2, is 
implemented in doing so. The results at t = t^ ina -| are compared with the 
exact solutions and with the solutions obtained by using the new grid system 
for all time, i.e., from t = 0 to t = fl - na -| without changing the grid 
system. 


3.3 Results and Discussion 

Various grid systems, as shown in Figs. 3.3 - 3.5, and initial 
conditions are implemented as described in Sec. 3.2. The results are 
compared with the exact solutions and the solutions obtained by using the 
second grid system for all time. For the case of exponentially stretching. 



the derivative contained in the Jacobian of transformation, J, x_, x , etc. 

S n 

can be computed analytically. Here, results using both analytically 
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Fig. 3.3. Grid systems for a circular geometry using 11 x 11 and 


21 x 21 grids. 


Fig. 3.4. Grid systems for a circular geometry using 21 x 21 and 


21 x 21 grids 


HI 

8 BM 

hH B is 

eMSsiiSSi 


Fig. 3.5. Grid systems for an elliptic geometry using 11 x 11 and 


21 x 21 grids. 







computed derivatives and numerically computed derivatives are shown. 

Before discussing any results, it should be noted that some error 
should be expected since Rams h aw' s technique assumes constant volumetric 
density of the conserved quantity in each grid cell. But in this study, the 
analytic initial conditions which vary within grid cells have been selected. 
However, the objective of this study is to determine whether the errors grow 
significantly as Ramshaw's idea is implemented. 

Tables 3.1 and 3.2 show the root mean square errors at final time, 
t = t (both tables are at 200 time steps). If errors are defined as the 
difference between the calculated solutions and the exact solutions, the 
root mean square errors can be obtained as 

I 2 ~ 

rms. error = / Terror 
yj HTxTT 

where m x n are the number of cells contained in the grid system where the 
errors are computed. 

Tables 3.1 and 3.2 indicate that the technique does not introduce 
significant errors to the solutions. The errors introduced by the technique 
can be seen clearly by looking at the top row of Table 3.2. In the case of 
constant mesh size with linear initial conditions, the finite volume calcu- 
lation should yield the exact solutions. Note that multiple grid calcu- 
lations even yield smaller errors than the single grid calculation in some 
cases. Thus, it can be concluded that the technique does not introduce 
significant errors to the solutions, especially when it is implemented on 
two grid systems which contain grid cells of comparative sizes. 
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Table 3.1 Root Mean Square Errors, 2D Case 1 


initial condition*, q ■ Acos(jx)sin(fy> 

exact solutions, q ■ Acos < I< x-at > > sin < £<y-bt ) > 


1 grid vs. 
£ grid 

m x n 
1 grid 

m x n 
£ grid 

rms. error 
1 grid 

rms. error 
2 grid 

rms. error 
multiple grid 

constant 
mash siza 

11x11 

£1 x£l 

.0110307 

. 0027987 

. 0029005 

expo anal. 

11x11 

£1x21 

. 028259 

. 007452 

. 007656 

k*1.5 num. 

11x11 

21x£l 

.02815 

.007414 

. 00762 

Fig 3.3a vs. 3.3b 

11x11 

21x21 

.0216039 

.009214 

. 0102846 

Fig 3.4a vs. 3.4b 

£1x21 

21x21 

. 002740 

.009214 

. 009668 

Fig 3.5a vs. 3.5b 

11x11 

21x21 

.0241553 

. 0099474 

. 0108937 
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1 grid v». m x n m x n rms. error rm«. trror rma. trror 

2 grid 1 Qrid 2 grid 1 grid 2 grid multiple grid 

conet ant 11x11 21x21 .0000001 .0000000 .0009590 

mesh size 

expo anal. 11x11 21x21 .0067252 .0017874 .0020591 

k«1.5 num. 11x11 21x21 .0064324 .0016984 .0019984 

Tig 3.3b vs. 11x11 21x21 .0025948 .0052553 .0075165 

Fig 3.3a & 3.4a 21x21 21x21 .0006518 .0052553 .0053213 

Fig 3.5a vs. 3.5b 11x11 21x21 .0052474 .0042585 .0043775 










4. THREE-DIMENSIONAL HYPERBOLIC EQUATION 


Consider the three-dimensional hyperbolic equation. 


q + aq + bq + cq = 0 0 <x <1, 0 <y<l, 0 <z<l 

t x y z 

with the initial condition, 


q( x, y, z, 0) = f ( x, y, z) 

The exact solution of this equation can be found as 

q(x, y, z, t) = f(x-at, y-bt, z-ct) 

As mentioned in Sec. 3, this problem is well posed if one specifies the 
boundary conditions at the boundaries where the dot product ( a,b,c) • ~n > 

0, where (a,b,c) = ae +be“ +ce and fT is the unit normal at the boundaries. 

This problem is solved on the multiple grid system using the finite- 
volume approach along with the three-stage Runge-Kutta integration in time. 
In Sec. 4.1, the formulation of the 3D hyperbolic equation on the curvilin- 
ear grid system is described. The brief discussion on the implementation on 
the multiple grid is given in Sec. 4.2. Section 4.3 discusses the results 
from various grid configurations and initial conditions. Note that the 3D 
hyperbolic equation is more suitable to test Ramshaw' s technique. This is 
because most real physical computations are performed on a three dimensional 
space and the governing equations are similar to the 3D hyperbolic 
equation. 
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4.1 Formulation of 3D Hyperbolic Equation on Curvilinear Grid 
Consider the three dimension! hyperbolic equation, 

i 

i 

q t + aq x + bq y + cq z = 0 (4.1) 

which can be written in the conservative form as [46] 

q t + f x + V h z ■ 0 <4 - 2) 

where 

f = aq, g = bq, h = cq 

and a, b, c are constants. 

This section shows the formulation of the transformed equation of Eq. 

4.2 on the computational coordinates (s-n-c). The advantages of trans- 
forming Eq. 4.2 into the computational coordinates were mentioned in Sec. 

3.1. Figure 4.1 illustrates the physical versus computational domain. 

| 

By using the chain rule, Eq. 4.2 is written as 

"t * ( Vx + Vy * Vz> + ( V>x + Vy + Vz* 

* ( Vx + Vy + Vz> = °- 

After some mathematics manipulations, one can write ! 

(I q) t + (U/ + 15 y g + ^5 z h ) ? + 0* x f + ln y g + ln^ 

J J J J J J J 
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(4.3) 


where 


+ (I5 x f + i ? y 9 + = 0 

J J J 



Similar to Eq. 3.6, one can see that 





Vc' I"z = - ( V C - W’ i?z 
J J 



Equation 4.3 can be written as 
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(4.4) 


where 









= 0 


q = 


q 


j 


f = l_(5 x f +5y g + 5 2 h) = (y n z 5 
g = lh x f * n y g * -> z h) - -( y f z ; 
h = + t^ + c z h) = (y^ 


z _ y_ ) f - (x z - z x )g + (x y - y X )h 
n c t) c t) 5 n C n £ 


h\ )f * ( y c • h\ )g - ( v< ■ v c )h 


Z c y ) f - (x r z - z X )g + (x y - y x )h 
C n ; n £ n £ n £ n 


J 

Equation 4.4 is the transformed equation which is solved on the computa- 
tional coordinate (£,n,c). 

4.1.1 Difference Equation 

Figure 4.2 illustrates a volume element on the computational coordi- 
nates. The difference equations are derived by applying the finite volume 
approach to Eq . 4.4 with the volume element as in Fig. 4.2. The details of 
the formulation are given below. 

Recall Eq. 4.4, 


or 
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Integrating over the volume element, it is expressed as 


Applying the divergence theorem to the right 
can be written as 


+h )dV 

c 

hand side, the above equation 




-(( t 

[J A C 


* + 




dA) 


Utilizing Fig. 4.2, the difference equation for the above equation can be 
written as 


ACAnAC'*' 


<t 


i + JL» J + I, k + 1 
2 2 2 


= -AnAc(f. + 1§ j + J. f k + 1 

yv 2 2 

T i , j + U k + V 
2 2 

-aca c ( g i + j + i t k + 

2 2 

'% + J., j, k + j, ) 

2 2 

-AfAnC^ + i $ j + i t k + 1 
^ 2 2 
-h i + 1 , j + 1 , k^ 


(4.5) 


where ’ U, j + 1, k + 1 = [( Vs ‘ Vc )f " ( Vc ' Vc )9 


+ { V ? ’ Vc )h] i, j + l, k +i 

2 2 
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1 , j + 1, k + 1_ “ — ^ x i , j + 1, 
2 


2 2 


k ‘ x i, j, k 


An 


+ X i , j + 1, k + 1 " X i , j, k + 1^ 
An 


= \ ( y 

i , j + 1 , k + 1 — i > j + 1, 

2 2 2 An 


k " y i, j, k 


* y i, j + 1, k + 1 " y i, j, k +1 
An 


i , j + I* k + 1 " - Z i, j + 1. k ~ Z i, j, k 
2 2 2 An 

* z i, j + 1, k + 1 " Z i, j, k +1 } 
An 


= 1 ( X . . , , - X . . 

i , j + I» k + I - i. h k + 1 i, J, k 

2 2 2 An 

+ x i, j + 1, k+l~ X i,J + l, k^ 


An 


= 1 ( V - y ) 

i » j + i» k + JL — y i , j» k + 1 i , J + 1> k 

2 2 2 A? 


+ y i , j + l t k + 1 " y i , j, k 
A? 


i , j + U k + l ' — ^ Z i , j, k + 1 Z i, j, k 

2 2 2 A c 

* 2 i, j + 1, k + 1 ~ Z i , j + 1, k } 
A? 
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q i + i, j, k + 2 = [ 

2 2 


( Vc ' Vc )f * ( V< - h\ ] * 


' ( Vc - Vc )h] i + 1, J, k.i 

2 2 


= l(x 


- X 




1 I 

2 

2 AC 




+ x . , . 

, - x . 



i + 1, j, k + 

1 1, J, 



AC 



= Ky. . . . . 

- y. . , 

i + L 

j, k + 1 

- i + 1, j, k 

i, J, k 

2 

2 

2 AC 



+ y i 

+ 1, j, k + 1 

y i, J, k +1 



A? 




= 1 ( z . 

- z . 

l + 1, 

j, k + 1 

- 1 + 1, J, k 

1, J, k 

2 

2 

2 AC 



+ z . 


z . 


i 

+ 1, j, k + 1 

1, j, k +1 



AC 




= l(x . . . , 

- X . 

1 + I, 

j, k + 1 

- 1 » J» k + 1 

i, J, k 

2 

2 

2 A£ 



+ X. 


X . 


1 

+ li j, k + 1 

1 + 1* j, 1 



AC 


i + l. 

j, k + 1 

= - y i, J, k + 1 

’ y i, J, k 

2 

2 

2 AC 



+ y - y ) 

i + 1, J, k + 1 + 1, j, k 


AC 
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i + h, j , k + 1_ = l^ z i , j, k + 1 ~ z i , j, k 

2 2 2 Ac 

+ z i + 1, j, k + 1 ~ z i + 1, j, k } 

AC 


/N 

h . 


f + i- J + b K * [<V„ - * e y„)f - < X E Z „ - V„>9 


( Vn - VnH *1, J + 1. k 

2 2 


= 1( x . , . , - x . 

i + _1 , j + 1., k - i +1, j, k i, j, k 
2.2 2 A? 

* X i + 1, j + 1, k ' X i, j + 1, k } 
AC 


yc 


= I ( y - y 

i + 1, j + 1, k - i + 1, j, k 1, j, k 

2 

2 2 AC 

+ y i + 1, j + 1, k " y i, j + 1, k^ 

AC 


1 + I, j + 1, k - 2 i + 1, j, k Z i, j, k 

2 2 2 AC 

+ Z i + 1, j + 1, k ' Z i, j + 1, k } 

AC 


*n 


= 1 ( X . . . , - X . . , 

i + i, j + i, k - i, j + 1, k 1, li k 

2 2 2 An 

* X i + 1, j + 1, k X i + 1, j, k^ 
An 
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i + 1, 
2 


j + U k ~^ y i, j + 1, k ~ y i, j, k 
2 2 An 

+ y i + 1, j + 1, k y i + 1, j, k^ 

An 


i + j., j + 1., k ~ - Z i, j + 1, k " Z i, j, k 

2 2 2 An 

+x. , ,,-x. 

i + 1, j + 1, k i + 1, j, k 

An 


f = aq, g = bq, h = cq 


q i, j + 1, k + 1 ~ 1 ^ q i + 1, j + 1, k + 1 + q i - 1, j + 1, k + 1 


) 




i + 1, j, k + 1 8 I (q i + 1, j + 1, k + 1 + q i + 1, j - 1, k + l 1 


q i + 1, j + 1, k = - (q i + 1, j + 1, k + 1 + q i + 1, j + 1, 


k - r 


q 


i + 1, j + k + 1^ 
2 2 2 


1 q 

7 


i + i» J + A> k + _1 
2 2 2 


1 

J 


= [x (y z - z y ) - x (y z -zy) 

i + 1, j + 1, k + 1 L 5 n c nV n ? C 5? 


+ x C ( Vn ' Vn n i *1, J ♦ 1, k+1 

2 2 2 
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i + l, j + 1, k + 1 - ^ x i + 1, j, k i, j, k 

- - - 2 2 

2 2 2 as 

+ X i + 1, j + 1, k " X i, j + 1, k^ 

AS 

+ “ ^ X i + 1, j, k + 1 X i, j, k +1^ 

2 A? 

* X i + 1, j + 1, k + 1 ~ X i, j + 1, k + 1 ^ 

AS 


*n 


= 1 [1( x . . - x . 

i + 1, j + 1, k + 1 - - i, J + 1, k i, J, k 

— — “22 

2 2 2 An 

* X i + 1, j + 1, k ~ X i + 1, j, k } 

An 

+ - (X i, j + 1, k + 1 ' X i, j, k +1 } 


2 
+ X 


An 

i + 1, j + 1, k + 1 " X i + 1, j, k + 1 )] 


An 


= 1 [ l(x. . . . - x . . , 

i + 1, J + 1, k + 1 - - i» J» k + 1 i, J. k 

— 2 2 
2 2 2 A? 

+ X i + 1, j, k + 1 X i + 1, j, k ^ 

AC 

+ “ ^ X i» J + 1, k+1 X i, j + 1, k ^ 

2 Ac 

* x i + 1, j + 1, k + 1 ~ x j + 1, j + 1, k )] 

Ac 


Similar expressions can be written for y , y , y , z , z , and z 

£ n ? t, n 5 

Eq. 4.5 can be written in the form: 


. Finally, 
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qt = f (q h ) 


The three-stage Runge-Kutta integration as described in Sec. 3.1.2 can be 
applied here. 

4.2 IMPLEMENTATION ON MULTIPLE GRID SYSTEM 

Figure 4.3 illustrates an example of multiple grid used in this study. 
Here, the grid system changes from one to another at the plane z = iy In 
order to solve the difference equations, as in section 4.1.1, on this grid 
system, the extra boundary conditions are needed at the interface, i.e., 
plane z = z^ where the two single grid system meet. These boundary condi- 
tions are not the physical boundary conditions. Many authors refer to them 
as interface conditions, for example see Berger [43]. Special treatments 
mush be given for these boundary conditions. Figure 4.4 shows the typical 
interface in this study. 

In Fig. 4.4, h Q denotes the flux h coming out from the first grid 
system and h^ denotes the flux h going into the second grid system. Here, 
h is needed for the first grid calculations and h is needed for the calcu- 
1 at ions on the second grid. They are the extra boundary conditions 
mentioned previously. In this study, h Q is obtained by the interpolations 
between the first and the second grid. The flux h n is obtained by trans- 
ferring h Q across the interface using the Ramshaw's technique described in 
Sec. 2. The calculations begin from initial time, t = 0, to some time, say 
t = t ir ^ na -|. The results from the multiple grid calculations are compared 
with the exact solutions and the solutions obtained from the single grid 
calculations for various grid system and initial conditions. 
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z 


Fig. 4.4. An interface where h 0 is obtained by 

interpolating h from the points indi- 
cated by dotted triangle and h n is ob- 
tained by Ramshaw's technique. 



4.3 RESULTS AND DISCUSSIONS 


Tables 4.1 and 4.2 show the results obtained from the calculations on 
various grid systems with two different initial conditions. If the error is 
defined as the difference between the numerical solutions and the exact 
solutions, the root mean square error can be found as 


rms. error 


m .li 


error 
R — 


where N = total number of grid cells. 

Both tables illustrate the rms. errors from the multiple grid calcu- 
lations as well as the calculations from each of single grid which compose 
the multiple grid. Table 4.1 gives the results with linear initial condi- 
tions, while the results with trigonomitrical initial conditions are given 
in Table 4.2. It can be seen from both tables that the multiple grid calcu- 
lation does not introduce significant errors to the over all solutions. In 
some cases, the multiple grid calculations even give better solutions than 
the calculations from the single grid. 


5. CONCLUSIONS 

The results from this study have shown that the technique works reason- 
ably well when the two grid systems contain grid cells of comparative sizes. 
This can be clearly seen when one recalls the assumptions made for the tech- 
nique, i.e., each grid cell contains uniform conserved quantity. Since, in 
general, this criterion cannot be met in the real physical calculations, 
some effort may need to be made to seek for the best possible types of grids 
in order to apply the technique. It should be pointed out that only simple 
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Table 4.1 Root Mean Square Errors, 3D Case 1. 


initial conditions, q ■ A+x+y+z 

exact solutions, q ® A+ (x-at ) + (y-bt ) + ( z-ct ) 



No. of pts. l*^vs. £** 


rms 

. error 

conf i gurat ion 

single grid 

mul. grid 



first 

second 



llxllxll vs. £1x21x11 

. 000000 

. 000000 

. 013660 


£lx£lxll vs. llxllxll 

. 000000 

. 000000 

. 000004 


17x17x11 vs. £1x21x11 

. 000000 

. 000000 

. 005720 

const, mesh 

£1x21x11 vs. 17x17x11 

. 000000 

. 000000 

. 003334 


19x19x11 vs. 21x21x11 

. 000000 

. 000000 

.002163 


21x21x11 vs. 19x19x11 

. 000000 

. 000000 

.001703 

exp. streching 

21x21x11 vs. 21x21x11 

. 000000 

. 002106 

. 002389 

w i th k=3. O 

21x21x11 vs. £1x21x11 

. 002106 

. 000000 

.007139 

( *=const . mesh > 
cy 1 i rider 

□-type vs. H-type 
£1x21x11 vs. 21x21x11 

. 000608 

.015882 

. 009620 


llxllxll vs. 21x21x11 

. 002842 

.012321 

. 012883 

But ler-wing 

17x17x11 vs. 21x21x11 

.001495 

.012321 

. 008672 


21x21x11 vs. 21x21x11 

. 001275 

.012321 

. 008379 
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Table 4.2 Root Mean Square Errors, 3D Case 2 


initial conditions, q * flcos (lx ) sin (fy ) sin <f z) 

exact solutions, q « Ocos ( I(x-at ) ) sin (I (y-bt ) ) sin ( \ (z-ct ) > 

tit 


conf i gurat ion 

No. of pts. l^vs. 2*^ 


rms 

. error 



single grid 

mul. grid 



| 

m 

second 



llxllxll vs. 21x21x11 

. 005560 

.001901 

.019420 


21x21x11 vs. llxllxll 

. 001901 

. 005560 

. 003786 

const, mesh 

17x17x11 vs. 21x21x11 

. 002494 

.001901 

. 006363 


21x21x11 vs. 17x17x11 

.001901 

. 002494 

.004219 


19x19x11 vs. 21x21x11 

. 002136 

.001901 

.003414 


21x21x11 vs. 19x19x11 

.001901 

.002136 

. 003070 

exp. streching 

21x21x11* vs. 21x21x11 

. 001901 

. 004424 

. 005578 

with k*3. 0 

21x21x11 vs. 21x21x1 1* 

. 004424 

.001901 

. 008227 

( *=const . mesh ) 
cyl inder 

O-type vs. H-type 
21x21x11 vs. 21x21x11 

. 001651 

.021766 

.011515 


llxllxll vs. 21x21x11 

. 007848 

.0211753 

. 017058 

But ler-wing 

17x17x11 vs. 21x21x11 

. 003297 

.0211753 

.012327 


21x21x11 vs. 21x21x11 

.002180 

.0211753 

.011819 
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configurations are chosen for this study, i.e., a single grid system can be 
generated around these configurations. For more complex bodies, the 
multiple grid approach may be the only way to attack the problem if the 
complexity of the geometries is to be maintained. 

This study demonstrates that use of the multiple grid does not intro- 
duce any significant error to the problems. The next step is to apply this 
technique to solve the equations of motion over an ideal aircraft config- 
uration such as a Butler-Wing configuration (Fig. 5.1). If the results from 
the Butler-Wing calculations are satisfied, the next step is to 
consider this technique along with the equations of motion over the real 
aircraft configuration. It is expected that results from this study will 
give some contributions to the field of Computational Fluid Dynamics. 
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